#define AnaND280HL2_cxx
#include "AnaND280HL2.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TFile.h>

#include "ND280HL2Utils.h"

void AnaND280HL2::DoAnalysis(const TString filelist, const Int_t ntarget, const TString tag)
{
  NeutrinoUtils::IniCC();

  Init(ND280HL2Utils::IniIO(filelist, "default"));

  printf("AnaND280HL2::DoAnalysis ntarget %d\n", ntarget);

  TList *lout=NeutrinoUtils::GetHistList();

  TTree *tout = NeutrinoUtils::GetTree();

  const Int_t nentries = fChain->GetEntriesFast();
  Int_t nsel=0;
  Int_t nfail = 0;
  for (Int_t jentry=0; jentry<nentries;jentry++) {
    if(jentry%10000==0){
      printf("%d/%d nsel %d nfail %d\n", jentry, nentries,nsel, nfail);
    }

    if(!fChain->GetEntry(jentry)){
      printf("AnaND280HL2::DoAnalysis GetEntry ends at jentry %d\n", jentry);
      break;
    }

    //--------------------------------------------
    NeutrinoUtils::FillCount(lout, 0);

    ND280HL2Utils::SaveInfo(nu_truereac, nu_parentpdg, nu_pdg, ntpctracks);

    /*
    const Int_t nhitcut = -999;
    if(NeutrinoUtils::fLeptonNhits[0]< nhitcut || NeutrinoUtils::fLeptonNhits[1] < nhitcut){
      nfail++;
      continue;
    }
    NeutrinoUtils::FillCount(lout, 2);
    */

    //=====================================================================
    //with NDHL2, only selected/reconstructed muon per event is save
    //simulated muon can be null (NC)
    //=====================================================================

    //---------------------------------------------------------------------------------------------------------
    //reconstruction
    //---------------------------------------------------------------------------------------------------------
    ND280HL2Utils::SetNeutrinoRec(vertex_pos);

    const Int_t nmurec = ND280HL2Utils::SetParticleRec(NeutrinoUtils::fNeutrinoRec, NeutrinoUtils::fMuonRec, NeutrinoTools::MuonMass(), selmu_dir, selmu_mom[0], vertex_pos, selmu_endpos, selmu_Chi2[0], selmu_NDOF[0], selmu_ntpcs, selmu_nfgds, selmu_tpc_nhits, selmu_tpc_charge);
    NeutrinoUtils::SetMuonNRec(nmurec);

    /*
    //allow multiple muon
    if(NeutrinoUtils::GetMuonNRec()>1){
      nfail++;
      continue;
    }
    NeutrinoUtils::FillCount(lout, 4);
    */

    /*
    //allow 0-reconstructed in order to obtain efficiency
    if(NeutrinoUtils::GetMuonNRec()==0){
      nfail++;
      continue;
    }
    NeutrinoUtils::FillCount(lout, 5);
    */

    ND280HL2Utils::SetNeutrinoGuess(nu_truedir);
    ND280HL2Utils::SetMuonGuess();

    NeutrinoUtils::SetRecGuessKinematics(kTRUE);

    //---------------------------------------------------------------------------------------------------------
    //simulation 
    //---------------------------------------------------------------------------------------------------------

    ND280HL2Utils::SetNeutrinoSim(nu_truedir, nu_trueE);

    //test understanding of coordinate system ->
    //const TVector3 * ND280HL2Utils::GetNuDirSim(const Float_t nu_NuParentDecPoint[], const Float_t vertex_truepos[], const Float_t nu_truedir[])
    //const TVector3 * testtmp = ND280HL2Utils::GetNuDirSim(nu_NuParentDecPoint, vertex_truepos, nu_truedir);
    //delete testtmp;
    //test<-

    const Int_t nsim = ND280HL2Utils::SetParticleSim(NeutrinoUtils::fNeutrinoSim, NeutrinoUtils::fMuonSim, NeutrinoTools::MuonMass(), truemu_truedir, truemu_truemom, particle);
    NeutrinoUtils::SetMuonNSim(nsim);

    NeutrinoUtils::SetSimKinematics();
    /*
    //with single CC muon
    //open to allow multi-(+/-)muon real case: decay into muon and numu
    if(NeutrinoUtils::GetMuonNSim()>1){
      nfail++;
      continue;
    }
    NeutrinoUtils::FillCount(lout, 6);
    */

    /*
    //allow NC
    if(NeutrinoUtils::GetMuonNSim()==0){
      nfail++;
      continue;
    } 
    NeutrinoUtils::FillCount(lout, 7);
    */

    //---------------------------------------------------------------------------------------------------------
    //---------------------------------------------------------------------------------------------------------
    NeutrinoUtils::FillCount(lout,11);

    tout->Fill();
    
    nsel++;

    if(ntarget>=0 && nsel>=ntarget){
      break;
    }
  }

  printf("nentries %d nsel %d nfail %d delta %d\n", nentries, nsel, nfail, nentries-nsel-nfail);

  TFile *fout=new TFile(Form("ND280HL2AnaOutput_tag%s_%d_%d.root", tag.Data(), nsel, nfail),"recreate");

  tout->Write();
  //for further handling in plot
  lout->Write("lout",TObject::kSingleKey);

  fout->Save();
  fout->Close();
  delete fout;
}
